\section{三维重建理论基础}\label{preliminaries}

本章将介绍三维重建相关的理论，首先从基础的数学几何模型入手对三维重建问题进行建模，接着介绍了传统三维重建算法的两个重要部分：运动恢复结构、多视图立体几何。最后介绍了端到端深度学习三维重建。

\subsection{三维重建基础理论}

本节主要讲解三维重建整个流程中比较重要的数学几何理论，包括相机模型与相机坐标系、对极几何与极线约束、单应性变换。


\subsubsection{相机模型与坐标系}\label{sec:axis}
本小节主要需要讨论的问题是：不同相机拍摄的三维空间点是如何投影到二维图像上的。这实际上是一个坐标系转换的问题， 在三维重建过程中，常讨论的坐标系包括四个：世界坐标系、相机坐标系
、图像坐标系、像素坐标系。

不同相机成像方式不同，但在相机坐标系与世界坐标系上是保持一致的，差异主要在于图像坐标系与像素坐标系。

\begin{enumerate}
    \item \textbf{针孔相机模型与坐标系}:
    
    对于针孔相机来说，四个坐标系的关系如下：

    \begin{itemize}[leftmargin=4em]
        \item \textbf{世界坐标系:} 以某一参考点构建的坐标系。所有相机都相对该参考点构成坐标系，如图\ref{fig:axis:c}。
        \item \textbf{相机坐标系:} 根据刚体运动原理将世界坐标系旋转平移到相机坐标系，如图\ref{fig:axis:b}。 
        \item \textbf{图像坐标系:} 以图像中心为原点，同时单位不再是像素而是其缩放，此时成像平面位于距离相机镜头一个焦距长度的位置，如图\ref{fig:axis:a}。
        \item \textbf{像素坐标系:} 以图像左上角为原点，横向为u轴，纵向为v轴的二维坐标系。这个坐标系下的图像可以理解为程序直接读取到的图像数据，即二维数组，如图\ref{fig:axis:a}。
    \end{itemize}
    \begin{figure}[H]
        \centering
        \subfigure[图像坐标系与像素坐标系]{
        \begin{minipage}[t]{0.25\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/axis1.png}
        \end{minipage}\label{fig:axis:a}
        }
        %
        \centering
        \subfigure[相机坐标系]{
        \begin{minipage}[t]{0.25\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/axis2.png}
        \end{minipage}\label{fig:axis:b}
        }
        \centering
        \subfigure[世界坐标系]{
        \begin{minipage}[t]{0.25\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/axis3.png}
        \end{minipage}\label{fig:axis:c}
        }
        \caption{相机坐标系}\label{fig:axis}
        \centering
    \end{figure}
    
    
    上述关系最终用公式表述为：
    
    \begin{equation}
        \begin{pmatrix}
            u\\v\\1
        \end{pmatrix}
        =\frac{1}{z_{c}}\begin{pmatrix}
            f_{x}&0&c_{x}\\0&f_{y}&c_{y}\\0&0&1
        \end{pmatrix}
        \begin{pmatrix}
            x_{c}\\y_{c}\\z_{c}
        \end{pmatrix}
        =\frac{1}{z_{c}}KP_{c}=p
    \end{equation}
    
    公式中 $f$表示焦距，$c$表示图像中心到图像边界的距离。$\frac{1}{z_{c}}$表示的是尺度，也是希望从图像中恢复出来的重要的量。上述公式描述了从世界坐标系下真实点到相机图像上的映射。
    
    
    \item \textbf{全景相机模型与坐标系}
    
    全景图的相机模型可以视作一个球形(如图\ref{fig:e2p})，则全景图的x轴与y轴可以转化为对应的经度$\lambda$与纬度$\varphi$：
    
    \begin{figure}[H]
        \centering
        \subfigure[相机图像坐标系]{
        \begin{minipage}[t]{0.4\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/sphere.png}
        \end{minipage}\label{fig:det}
        }
        %
        \centering
        \subfigure[像素坐标系]{
        \begin{minipage}[t]{0.44\linewidth}
        \centering
        \includegraphics[width=0.8\linewidth]{figures/perspective.png}
        \end{minipage}\label{fig:seg}
        }
        \caption{全景相机模型}\label{fig:e2p}
        \centering
    \end{figure}
    
    \begin{equation}
        \begin{split}
            \lambda &=\frac{x}{W}2\pi \\
            \varPhi &=\frac{y}{H}\pi      
        \end{split}
    \end{equation}
    
    其中$W$为图片宽度，$H$为图片高度，$\varPhi$为纬度，$\lambda$为经度。同时球体的$xyz$坐标(相机坐标系)可以与经纬度坐标进行转化：
    
    \begin{equation}
        \begin{split}
           x&=r\cdot \cos (\varphi)\sin(\lambda)\\
           y&=r\cdot \sin (\varphi)\\
           z&=r\cdot \cos (\varphi)\cos(\lambda)
        \end{split}
    \end{equation}
    
    这里$r$表示球坐标的半径，也即深度。
    
    \begin{figure}[htbp]
        \centering
        \includegraphics[width=0.7\textwidth]{figures/FOV.png}
        \caption{视角}
        \label{fig:FOV}
    \end{figure}
    
    最后给出视角的定义：视角是一个相机可以成像的最大范围，它与两个因素有关：镜头的焦距和传感器的大小(如图\ref{fig:FOV})。其关系如下：
    \begin{equation}
        AFOV=2\tan^{-1}(\frac{h}{2F})
    \end{equation}
    
    这里$AFOV$即角度制的视角大小，$F$为相机焦距，在全景图中可以为转换后的焦距。$h$为图像高度。
    
    \textbf{全景相机转针孔相机模型}\label{sec:e2p}
    在给出$FOV$、$W$、$H$的条件下，可以推算出Focal、$c_{x}$、$c_{y}$，则得到待求针孔相机模型的内参矩阵$K$。则对于待求图像的任意一点($u$,$v$,1)，可以转换到相机坐标系($x$,$y$,$z$)：
    
    \begin{equation}
        \begin{pmatrix}
            x\\y\\z
        \end{pmatrix}=\begin{pmatrix}
            u\\v\\1
        \end{pmatrix}\begin{pmatrix}
            f_{x}&0&c_{x}\\0&f_{y}&c_{y}\\0&0&1
        \end{pmatrix}^{-1}
    \end{equation}
    这里 $u$、$v$为图像坐标系中的坐标。
    
    接着根据相机坐标系与经纬坐标系的关系、经纬坐标系与全景图像素坐标系的关系，可以在透视图像素坐标系与全景图像素坐标系之间建立等式。即可完成该转化。
    
    
\end{enumerate}

\subsubsection{单应性变换}

单应性变换是一种二维到二维的投影，又名射影变换，也可以理解为从像素坐标系到像素坐标系的变换。
该变换可以用一个非奇异的矩阵H表示。假设相机A与相机A'之间的相对位姿为[$R$|$t$]，则公共点x其在相机坐标系下存在关系：

\begin{equation}
    x^{'}=Rx+t
\end{equation}

带入像素坐标系与相机坐标系的转换关系后，可以得到单应性变换矩阵H为：
\begin{equation}
    H=K(R+\frac{1}{d}tn^{T})K^{-1}
\end{equation}

这里$n$为法向量。$R$、$t$为相机位姿，$K$为内参矩阵，$d$为深度。

在代数中，图像上一点p=(x,y,1)通过单应性矩阵转换到另一点p'=(x',y',1)，可以得到$p'=Hp^{T}$,将矩阵H展开为九维向量可以得到：
\begin{equation}
    \begin{pmatrix}
        -x&-y&-1&0&0&0&xx'&yx'&x'\\
        0&0&0&-x&-y&-1&xy'&yy'&y'
    \end{pmatrix}h=0
\end{equation}

规定左侧矩阵为$A$，由于单应性矩阵自由度为8，则只需要找到4个点对构成8维矩阵即可解得$H$，如果点对较多可以考虑通过最小二乘，对$A$矩阵进行SGD分解即可求解。

\subsubsection{对极几何}\label{sec:duijijihe}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/duiji.png}
    \caption{对极几何}
    \label{fig:duiji}
\end{figure}

对极几何（Epipolar Geometry）是对立体视觉建模的一种约束，通过此建模方法使得立体匹 配或深度估计问题存在一个最优解。如图\ref*{fig:duiji}所示，
P为空间中的物点，可以被相机$C_{l}$与$C_{r}$看到。$C_{l}$与$C_{r}$为相机中心，$I_{l}$与$I_{r}$为图像平面，$p_{l}$与$p_{r}$为物点p在两图像中的投影点。
相机中心连线为基线，左摄像机在右图中的投影点与右摄像机在左图的投影点为极点。相机原点、物点构成的平面为极平面。极平面与两个相平面交线为极线$l_{pr}$、$l_{pl}$。而极线约束即：
给定图像上的一个特征，它在另一幅图像上的匹配视图一定在对应的极线上，即已知$p_{r}$,则它对应在右图的匹配点一定在极线$l_{pr}$上反之亦然。

极线约束给出了对应点的重要的约束条件，它将对应点匹配从整幅图像中查找压缩到一条线上查找，
大大减小了搜索范围，对对应点的匹配起指导作用。

假设第一帧的相机坐标系下，点P的空间位置为$P=[X,Y,Z]^{T}$。则P在第一帧的像素坐标为$p_{1}$，第二帧为$p_{2}$。则有关系：

\begin{equation}
    s_{1}p_{1}=KP,s_{2}p_{2}=K(RP+t)
\end{equation}

$K$为内参矩阵，$R$、$t$为相机的相对运动。由于最终构建的是齐次方程，所以可以省略$s_{1}s_{2}$，则像素点在归一化平面上的坐标为：

\begin{equation}
\begin{split}
    x_{1}&=K^{-1}p_{1}\\
    x_{2}&=K^{-1}p_{2}\\
    x_{2}&=Rx_{1}+t
\end{split}
\end{equation}
等式同时对t做外积，同时左乘$x^{T}_{2}$，可以得到:

\begin{equation}
    \begin{split}
        &t\times x_{2}=t\times Rx_{1}\\
        &\Leftrightarrow x^{T}_{2}t\times x_{1}=x^{T}_{2}t\times Rx_{1}\\
        &\Leftrightarrow x^{T}_{2}t\times Rx_{1}=0\\
        &\Leftrightarrow p_{2}^{T}K^{-T}t\times RK_{-1}p_{1}=0
    \end{split}
\end{equation}

这便是对极约束，几何含义是两个相机原点与物点三者共面。将最后两行公式的中间部分分别提炼出来，记作基础矩阵（Fundamental Matrix）F与本质矩阵（Essential Matrix）E。

\begin{equation}
    \begin{split}
    E&=t\times R\\
    F&=K^{-T}EK_{-1}
\end{split}
\end{equation}

\subsection{运动恢复结构}

三维重建输入了一系列无序的图片，重建空间的第一步是恢复每张图片的对应的空间视角。空间视角可以通过[$R$,$t$]表示，所以首先要求解的就是图片的[$R$,$t$]，把[$R$,$t$]
称为相机外参，与相机的焦距等参数构成的内参矩阵相对应。
运动恢复结构（Structure From Motion，简称SFM）的目标即对[$R$,$t$]进行解算，对空间坐标点进行定位，大体上可以分为三部分：特征点提取、特征点匹配、位姿解算与优化。
\subsubsection{特征点提取}
特征点提取是对图像中信息丰富、有价值的像素点进行提取与表示。在诸多特征点算法中
，以尺缩不变性变化（SIFT\cite[]{ng2003sift}）、加速稳健特征（SURF\cite[]{bay2006surf}）、特征点检测算子（ORB\cite[]{rublee2011orb}）为代表的特征描述算法取得了较好的效果。这些算法各有优缺点，

\begin{itemize}[leftmargin=4em]
    \item SIFT算法在处理尺度、旋转变化大的场景具有优势，对光照、噪声也有很好的鲁棒性，但是其运算量较大，实时性不足；
    \item SURF的计算量小于SIFT，同样具有尺度不变、旋转不变的特性，但实时性仍然不足；
    \item ORB运算速度快，具有旋转不变性，但是却不具备尺度不变性。
\end{itemize}

这些基于局部像素点的特征提取算法在弱纹理区域内并没有取得较好的效果\cite[]{张易2018弱纹理环境下视觉里程计优化算法研究}\cite[]{张庆鹏2021室内场景下弱纹理物体三维重建算法的研究}，
这是由于弱纹理区域在局部表现趋于一致，因此无法提取到较好的特征点。

另一方面，基于深度网络模型的特征点提取算法逐渐兴起，诸如Super Point\cite[]{detone2018superpoint}、LoFTR\cite[]{sun2021loftr}等。
SuperPoint采用自监督的方式，首先通过构建理想图形数据集，进行特征点提取，再利用单应性变换对输入的图片进行增广，最后通过理想图形训练的模型对输入图片及其变换进行预测，将结果作为真值进行训练。
SuperPoint的效果达到了SIFT等非深度网络模型特征点算子的效果，在某些指标上甚至超过了他们。
LoFTR则将NLP中的Transformer模型，引入注意力机制，使得特征点的提取不光关注于局部特征，同时关注特征点在全局的表现。这更符合人类对图像的认知规则。这些算法在弱纹理区域的特征点提取优于SIFT等传统算法。

由于本文选取的参考算法Colmap的特征提取为SIFT，故在这里详细介绍SIFT的算法原理。
SIFT算法主要特点是对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性。其主要流程为：从图像中构建差分高斯金字塔、
特征点定位、通过梯度确定方向、通过关键点周围梯度提取描述子。

\textbf{差分高斯核（DOG）：}DOG是一种常用的边缘和特征点检测方法。 将图像分别输入到$\sigma$不同的高斯核中，再将两个输出结果作差，DOG可以表示为：

\begin{equation}
    DoG\triangleq G_{\sigma_{1}}-G_{\sigma_{2}}=\frac{1}{\sqrt[]{2\pi}}(\frac{1}{\sigma_{1}}e^{-(x^{2}+y^{2})/2\sigma_{1}^{2}}-\frac{1}{\sigma_{2}}e^{-(x^{2}+y^{2})/2\sigma_{2}^{2}})
\end{equation}

不同尺度的高斯核表示了不同的模糊情况，二维高斯函数如图\ref{fig:guassian}所示:

\begin{figure}[htbp]
    \centering
    \subfigure[$\sigma$为2时的高斯函数]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[width=0.8\linewidth]{figures/sigma2_gs.png}
    \end{minipage}
    }
    %
    \centering
    \subfigure[$\sigma$为4时的高斯函数]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[width=0.8\linewidth]{figures/sigma4_gs.png}
    \end{minipage}\label{fig:del2}
    }
    \caption{二维高斯函数}\label{fig:guassian}
    \centering
\end{figure}

可以看到上述两个高斯核的形状不同，$\sigma$小的高斯核在中心更加集中，在边缘衰减较快，卷积之后模糊效果较弱。使用不同的高斯核可以达到不同的模糊效果以模拟人眼在不同尺度下的观测效果。

通过对像素点在DOG空间上同一张图片的邻域以及相邻尺度图片的邻域进行极值搜索可以得到特征点的位置。

\textbf{图像金字塔：} 金字塔即每层分辨率减少固定比率（图像大小在变，比如说原图像是640*320，变成320*160，再变成160*80），具有生成较快，且占用存储空间少（单纯降采样）的优点，而差分高斯核在进行高斯二维卷积的时候会产生冗余信息，图像金字塔则刚好可以弥补这一点。
SIFT将二者进行结合，即将图片进行图像金字塔，每层金字塔进行特定数量的高斯模糊，再进行差分。最终得到DOG金字塔。从DOG金字塔检测到的关键点通过插值可以获得精准坐标。

\textbf{边缘效应消除：} DOG算子有较强的边缘相应，需要取出不稳定的边缘点。主要方法是计算特征点Heesian矩阵(如公式\ref*{equ:hessian})的特征值。Hessian矩阵实际上是二维变量的二阶导数，用来表示各个方向上的梯度变化。其几何意义就是，其两个特征值越大，该点所在局部区域凸性越强，变化就越大。而边缘点则是一个特征值很小，而另一个特征值较大。
Hessian矩阵的特征值代表了x和y方向上的梯度。较大特征值的特征向量是垂直于直线的，较小的则是沿着直线方向的。而边缘恰好是这种情况，所以如果二者比率大于特定值则为边缘。

\begin{equation}
    H(f)=\begin{pmatrix}
        \frac{\partial^{2}f}{\partial x^{2}}& \frac{\partial^{2}f}{\partial x \partial y} \\
        \frac{\partial^{2}f}{\partial x \partial y}& \frac{\partial^{2}f}{ \partial y^{2}} \\ 
    \end{pmatrix}
    \label{equ:hessian}
\end{equation}

\textbf{特征点描述子：} SIFT的提取的特征（关键点）需要保持尺度不变性，因此最终的特征点需要携带尺度信息，包括缩放、方向等信息，主要由它周围的像素来贡献。如果用一个向量来存储这些信息，此向量称为关键点描述子。
生成关键点的描述子之前，要先求出该关键点的主方向。为了让关键点对方向具有不变性，在选取16x16矩形区域的时候，将矩形的方向旋转到与关键点主方向一致。为了使描述符具有旋转不变性，需要利用图像的局部特征为给每一个关键点分配一个基准方向。使用图像梯度的方法求取局部结构的稳定方向。

SIFT定义关键点主方向为：以关键点为中心的周围像素所贡献的主方向。可采用梯度直方图统计法，统计以关键点为原点，一定区域内的图像像素点对关键点方向生成所作的贡献，贡献最大的那个方向即为关键点主方向。

具体方法如下：
\begin{enumerate}
    \item 对于在DOG金字塔中检测出的关键点，采集其所在高斯金字塔图像3$\sigma$邻域窗口内像素的梯度和方向分布特征；
    \item 在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0\~360度的方向范围分为36个区间，构建直方图。直方图的峰值方向代表了关键点的主方向(为简化，图中只画了八个方向的直方图)；
    \item 方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。这样在求解该关键点描述符时，取得16*16的矩形框以主方向为准；
    \item 为了增强匹配的鲁棒性，保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值。至此，将检测出的含有位置、尺度和方向的关键点即是该图像的SIFT特征点。
\end{enumerate}

经过上述过程，特征点的所有量$（x,y,\sigma,\theta）$都已经已经求得，其中位置（x,y）、
尺度$\sigma$是在DOG金字塔中求得，而特征点方向$\theta$是通过特征点邻域直方图求得。
接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等等。
这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。具体流程如下：

\begin{enumerate}
    \item 在关键点所在图像上，划出以关键点为中心的16x16的矩形图像；
    \item 将16x16矩形图像划分为16小格，每小格为4x4，并计算每个像素的梯度和幅度（即像素值变化的方向及大小）；
    \item 对每个小格进行统计，统计8个方向的幅度，形成梯度直方图；
    \item 将16小格的幅度直方图连接起来，用向量表示，即为关键点描述子，共有128（8x16）维；
    \item 去除光照变化的影响，需对上述生成的特征向量进行归一化处理。
\end{enumerate}

这样即可形成128维的描述子。

\subsubsection{特征点匹配}

在特征点匹配方面，FLANN\cite[]{muja2009flann}、2NN、brute-force算法结合RANSAC算法取得了很好的鲁棒性。然而这些匹配算法匹配时并没有考虑特征的结构相似性以及外观相似性，在室内场景表现无法让人满意。
SuperGlue\cite[]{sarlin2020superglue}引入了Transformer模型与注意力机制，提高了效果，在室内数据集的表现超过了传统的特征点匹配算法。

本文选择参考的Colmap的匹配算法包含多种匹配方式，在这里选择主要部分进行介绍：

\textbf{匹配准则：}Colmap匹配特征点选用的是余弦相似度。余弦相似度用向量空间中两个向量夹角的余弦值作为衡量两个个体间差异的大小。余弦值越接近1，就表明夹角越接近0度，也就是两个向量越相似，这就叫"余弦相似性"。其计算方式如公式\ref*{equ:cos}

\begin{equation}\label{equ:cos}
    cos(\theta)=\frac{\sum_{i=1}^{n} (X_{i}\cdot y_{i})}{\sqrt[]{\sum_{i=1}^{n}(x_{i})^{2}} \cdot \sqrt[]{\sum_{i=1}^{n}(y_{i})^{2}}}
\end{equation}

\textbf{暴力匹配：} 暴力匹配即选择图片两两进行匹配，理论上是效果最好，但最耗时的。在实际操作时却因为图片的相似、重复结构导致效果不一定是最好。


\textbf{序列匹配：} 序列匹配要求输入图像是按照顺序的，这也比较符合实际使用情况。因为一般去采集图像时，都是按时序去采集，甚至直接从视频提取帧。
这种匹配方式会选择相邻n张图像匹配，往往效果会比较好。但是也会有其他的问题，比如在拍摄开始和拍摄结束，都拍摄的同一个物体，但由于误差累计，可能会重建出两个一模一样的（只是会有些偏差）。

\textbf{空间匹配：}—这种方式是利用了地理位置信息，也就是每张图像必须自带位置信息，例如gps。这种方式需要设置在多大范围内进行匹配，因此需要事先知道大概重建场景大小，以便选择合适参数。


\subsubsection{位姿解算与参数优化}\label{sec:sfm}



在获取匹配的点对之后，SFM算法对相机参数与空间建构进行重构。而重构方案又包含了增量式与全局式两种。增量式与全局式在于计算对象选取的不同，前者是递增地进行计算，而后者是
对整体图像进行计算，二者原理上相通，本节主要介绍共同的原理。

\textbf{初始化：} 首先需要选择优质的图片作为初始的图像对。选择的规则如下：

\begin{itemize}
    \item focal length 存在；
    \item 匹配点数从大到小排序，数量靠前的；
    \item 第二张候选影 像的匹配数目大于设定的阈值。
\end{itemize}

接下来希望进行从图像到图像的相对定位，计算位姿。参考\ref*{sec:duijijihe}节，位姿估计问题实际由两步组成：

\begin{enumerate}
    \item 根据匹配点的像素位置，求出E或者F；
    \item 根据E或者F求出R，t。
\end{enumerate}

限于篇幅，选择形式更为简单的E矩阵进行求解的介绍。

E矩阵含有9个未知数，在计算机多视图几何中\cite[]{hartley2003multiple}，可以证明的是该矩阵的特征值为[$\sigma$,$\sigma$,0],且该矩阵满足对极约束，
故该矩阵一共存在五个自由度。理论上
五对点即可完成求解
，但在实际操作中往往只考虑对极约束减少到自由度，
故选择八对点进行求解，即八点法。考虑一对点$x_{1}$=$[u_{1},v_{1},1]$,$x_{2}$=$[u_{2},v_{2},1]$。根据
对极约束有：

\begin{equation}
    \begin{pmatrix}
        u_{1}&v_{1}&1
    \end{pmatrix}=\begin{pmatrix}
        e_{1}&e_{2}&e_{3}\\
        e_{4}&e_{5}&e_{6}\\
        e_{7}&e_{8}&e_{9}\\
    \end{pmatrix}    \begin{pmatrix}
        u_{1}&v_{1}&1
    \end{pmatrix}^{T}
\end{equation}

可以将E展开，写为向量的形式,相应的$x_{1}$$x_{2}$也可以写成向量的形式：

\begin{equation}
    \begin{split}
        \overrightarrow{e}&=\begin{pmatrix}
            e_{1}&e_{2}&e_{3}&
            e_{4}&e_{5}&e_{6}&
            e_{7}&e_{8}&e_{9}&
        \end{pmatrix}^{T}\\
        \overrightarrow{vec_{1}}&=\begin{pmatrix}
            u_{1}^{1}u_{2}^{1}&u_{1}^{1}v_{2}^{1}&u_{1}^{1}&
            v_{1}^{1}u_{2}^{1}&v_{1}^{1}v_{2}^{1}&v_{1}^{1}&
            u_{2}^{1}&v_{2}^{1}&1&
        \end{pmatrix}^{T}\\
        \overrightarrow{vec_{1}}\cdot \overrightarrow{e}&=0
    \end{split}
\end{equation}

选取八对点即可构成一个线性方程组，如果该方程组满秩则可以求解。实际操作中，往往会选取多对点使用SVD进行最小二乘，求取优化解。
为了保证鲁棒性，会使用RANSC算法与之结合。

\textbf{后方交会（PnP）：} 经过初始化之后实际上获得了第一组图像的R、t以及对应特征点的空间坐标。而后续加入的图片如果继续使用2D-2D的算法
虽然仍然完成解算，但是却存在着需要点数较多、初始化、纯旋转和尺度的问题。因此希望可以利用已经获得的3D点的信息对相机的运动进行估计。这便是
3D-2D方法，而解决方法中比较经典的是构建Perspective-n-Point(PnP)问题（如图\ref{fig:pnp}）。最少只需要三对点即可完成解算。

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{figures/pnp.png}
    \caption{PnP}
    \label{fig:pnp}
\end{figure} 

\textbf{三角化：} 三角化的目的是获得图像中特征点的空间坐标，其原理如图\ref*{fig:tri}所示。

\begin{figure}[htbp]
    \centering
    \subfigure[]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[width=0.8\linewidth]{figures/tri1.png}
    \end{minipage}
    }
    %
    \centering
    \subfigure[]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[width=0.8\linewidth]{figures/tri2.png}
    \end{minipage}
    }
    \caption{三角化恢复空间坐标}\label{fig:tri}
    \centering
\end{figure}


\textbf{光束平差法（Bundle Adjustment\cite{triggs1999bundle}）：}中文名称是光束法平差,经典的BA目的是优化相机的pose和landmark,其在
SfM和SLAM 领域中扮演者重要角色。所谓的 BundleAdjustment，是指从视觉重建中提炼出最优的3D模型和相机参数（内参数和外参数）从每一个特征点反射出来的几何光束，在把相机状态和特征点空间位置做出最优调整之后最后收束到相机光心的这个过程。

假设$z_{i,j}$为相机在视角$v_{i}$处观察到的三维路标点$p_{j}$,从三维点到二维点的投影可以由\ref{sec:axis}节处推导，将其定义为$f$，则所有观测点与投影点的误差可以写为：

\begin{equation}
    \frac{1}{2}\sum_{i=1}^{m} \sum_{j=1}^{n}  \Vert e_{i,j} \Vert^{2}=\frac{1}{2}\sum_{i=1}^{m} \sum_{j=1}^{n} \Vert z_{i,j}-h(v_{i},p_{j})\Vert^{2}
\end{equation}

对上述最小二乘进求解则会对v与p同时做出调整。具体求解方法则涉及到非线性优化问题。主要求解方法包括：Gauss-Newton、Levenberg-Marquadt法及其衍生方法。需要解决的是R、t的求导问题。可以转化到李代数上从而完成求导运算。在这方面有很多优秀的开源库可以简化运算，包括ceres、g2o等。
\subsection{多视图立体几何}\label{sec:mvs}

多视图立体几何（MultiViewStereo）也即MVS是希望在给出稀疏点云与相机位姿的条件下将点云稠密化以进一步完成三维重建。
通常通过计算每一帧图片的深度图与法向量图，结合原始图像进行深度图融合得到稠密点云。
而根据深度图计算方式的不同，可以分出两个分支即Plane-sweeping与Patchmatch。
这两种算法都是希望获得图像的像素级深度图，再将其与图像进行融合得到稠密点云。
诸多开源库实现了这两种思路，包括Colmap\cite[]{schoenberger2016mvs}、openmvs\cite[]{openmvs2020}等。
围绕着这两个思路也延伸出了一系列端到端的深度学习算法，
包括基于Plane-sweeping的MVSNet系列，以MVSNet\cite[]{yao2018mvsnet}最为代表性，以及基于Patchmatch的PatchMatchNet\cite[]{wang2021patchmatchnet}。
因此Plane-sweeping与Patchmatch这两个算法极具参考意义，故在这里进行展开讲述。

在进行深度计算之前，需要在无序的图片中为每一个图片找到其关联性最好，最适合做深度计算的邻域帧。邻域帧的选择策略可以建模为一个马尔科夫随机场的Labeling问题，
由于与本文的研究目标关系不大在这里不做展开。

\subsubsection{Plane-sweeping}

如图\ref{PlaneSweeping},Plane Sweeping算法将深度范围划分为一个个平面，深度范围内可以由很多方法获得。如果平面足够密集，空间分割足够细腻，
那么，空间物体表面上的一点M一定位于众多平行平面中的其中一个平
面上。判断曲面上一点在哪个平面上的方法则是所有可以看到曲面上的M的相机看到的是同一个颜
色，也就是物体在点M本来的颜色。

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.5\textwidth]{figures/planesweeping.png}
    \caption{PlaneSweeping算法示意图}
    \label{PlaneSweeping}
\end{figure}



\textbf{深度图初始化：} 首先将SFM得到的稀疏点云进行Delaunay三角化，如图\ref{fig:del}。Delaunay三角具有如下性质：

\begin{itemize}[leftmargin=4em]
    \item 网格中的最小角最大化；
    \item 任意三角形的外接圆内不包含三角形以外的顶点；
    \item 最大化三角面片的内切圆的平均值。
\end{itemize}

\begin{figure}[htbp]
    \centering
    \subfigure[计算外接圆]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[width=0.8\linewidth]{figures/Delaunay1.png}
    \end{minipage}\label{fig:del1}
    }
    %
    \centering
    \subfigure[最终效果]{
    \begin{minipage}[t]{0.44\linewidth}
    \centering
    \includegraphics[width=0.4\linewidth]{figures/Del2.png}
    \end{minipage}\label{fig:del2}
    }
    \caption{Delaunay三角化}\label{fig:del}
    \centering
\end{figure}

随后按照三角形的三个点对深度图进行插值得到完整的初始化深度图，如图\ref*{fig:depthinit1}。

\begin{figure}[H]
    \centering
    \subfigure[图像三角化]{
    \begin{minipage}[t]{0.4\linewidth}
    \includegraphics[height=0.2\textheight]{figures/depthinit2.png}
    \centering
    \end{minipage}
    }
    \centering
    \subfigure[深度图插值]{
    \begin{minipage}[t]{0.4\linewidth}
    \centering
    \includegraphics[height=0.2\textheight]{figures/depthinit1.png}
    \end{minipage}
    }
    \caption{深度图初始化结果}\label{fig:depthinit1}
    \centering
\end{figure}

\textbf{匹配代价计算：}接下来需要计算每个像素点在不同深度层的匹配代价。匹配代价用于判断像素点是否在某深度平面上。常用的匹配代价包括NCC，Census等。
Census具有实现简单、速度较快的优势。Census基于像素邻域内的灰度差异，将像素灰度转换为比特串，即为Census值。在进行匹配时只需要计算Census的汉明距离即进行异或运算即可。

在参考图的一像素坐标$x$在当前深度下可以根据单应性矩阵H得到在右图的匹配点$x^{'}$。如果当前深度在真实曲面上，则$x$与$x^{'}$对应的像素点的表现是相似的，这种相似性可以通过匹配代价进行表达。
一张图片的匹配代价是x、y、d三维坐标下的函数，如图\ref{cost_volume}。

\begin{figure}[H]
    \centering
    \includegraphics[width=0.3\textwidth]{figures/cost_volume.png}
    \caption{匹配代价体}
    \label{cost_volume}
\end{figure}

\textbf{代价聚合：} 由于代价计算步骤只考虑了局部的相关性,对噪声非常敏感,无法直接用来计算最优视差,需要通过构建能量函数进行代价聚合对整体代价进行平滑。
通过该步骤,使聚合后的代价值能够更准确的反应像素之间的相关性。代价聚合最经典的算法是SGM，在速度和效果之间达到了平滑。SGM将二维的能量函数进行了一维近似，提高了效率。

在进行代价聚合之后每个像素直接选取聚合代价最小的视差值即是所求值。通常在得到视差值之后还需要进行左右视差图视察一致性的检查以对结果进行优化。

\subsubsection{PatchMatchStereo}\label{sec:mvs_pms}

PatchMatchStereo(PMS)算法的初衷是希望为所有像素点找到一个动态的视差平面，即视差平面根据物体表面而改变。因此为每一个像素都设置了一个视差与法向量，F=[(x,y,d),n]。每一个F都可以构成视差平面。
对于像素点p，其视差平面方程为：

\begin{equation}
    d_{p}=ap_{x}+bp_{y}+c
\end{equation}

其中a、b、c为视差平面的三个参数，$d_{p}$为视差，可以利用x、y与法向量n求得:
\begin{equation}
    \begin{split}
        a&=-\frac{n_{x}}{n_{z}}\\
        b&=-\frac{n_{y}}{n_{z}}\\
        c&=\frac{n_{x}x+n_{y}y+n_{z}d}{n_{z}}
    \end{split}
\end{equation}

所以PMS的目标就是找到每个像素的最优平面参数，即对每个像素找出聚合代价最小的那个平面。

\begin{equation}\label{equ:dispequ}
    f_{p}=arg \min_{f\in F}m(p,f)
\end{equation}

这里F是无边界的平面集合，m(p,f)为像素p当视差平面为f时的聚合代价值。计算公式为：

\begin{equation}
    m(p,f)=\sum_{q\in W_{p}}w(p,q)\cdot \rho(q,q-(aq_{x}+bq_{y}+c))
\end{equation}

这里$W_{p}$是一个以p为中心的方形窗口，w(p,q)是自适应的权值，为了解决边缘位置的错误匹配值问题（参考图\ref{fig:w}）。其原理是判断像素p和q在同一个平面的可能性。计算公式为：

\begin{equation}
    w(p,q)=e^{-\frac{\left\lVert I_{p}-I_{q} \right\rVert }{\gamma }}
\end{equation}
\begin{figure}[H]
    \centering
    \includegraphics[height=0.25\textheight]{figures/edge-flattening.png}
    \caption{权重w用来解决边缘误匹配问题}
    \label{fig:w}
\end{figure}


$\rho$ 用来衡量两个像素的不相似性。假设左视图像素q的视差平面参数为(a,b,c)，则它的视差值如式\ref*{equ:dispequ},在右视图上的同一点位$q'=q-d_{p}$,则二者的不相似性计算公式为：

\begin{equation}
    \rho(q,q')=(1-\alpha)\cdot min(\left\lVert I_{q}-I_{q'} \right\rVert,\tau_{col})+\alpha \cdot min(\left\lVert \Delta I_{q}-\Delta I_{q'} \right\rVert,\tau_{grad})
\end{equation}

公式中左项为像素强度之差，右项为梯度之差。其余参数为自定义参数，用以提高鲁棒性。

接下来介绍PatchMatchStereo流程：

\begin{enumerate}[leftmargin=4em]
    \item \textbf{随机初始：} 根据sfm的深度范围对每一个像素随机生成深度/视差和平面法向量；
    \item \textbf{迭代传播：} 迭代传播是PMS的核心所在，基本思想是把随机的所有视差平面中的少数正确的视差平面传播至同一视差平面内的其他像素。
    算法执行多次迭代，每次迭代执行四个步骤：
        \subitem （1）空间传播（Spatial Propagation）空间传播背后的思想是空间上相邻的像素极可能有相似的平面。基于此思想，PMS假设像素点p当前的平面为$f_{p}$，则检查p邻域内的像素q的平面$f_{q}$是否更适合p，即检查$m(p,f_{q})<m(p,f_{p})$,是否成立，
        若成立，则把平面$f_{q}$作为像素p的新平面。在偶数次迭代，q是p的左边和上边像素；在奇数次迭代，q是p的右边和下边像素。 
        偶数次迭代，从图像左上角像素沿行方向传播至右下角像素；奇数次迭代则和偶数次迭代传播顺序相反，即从右下角像素沿行方向传播至左上角像素。
        \subitem （2）视图传播（View Propagation）类似于一致性检查，PMS认为同名像对会有相似的视差平面。假设像素点p当前的平面为$f_{p}$，则在右视图上找到同名点为p的像素p′，
        $f_{p'}$为p′的视差平面转换至左视图坐标系后的平面，若则把平面$m(p,f_{p'})$<$m(p,f_{p})$，则把$f_{p'}$作为像素p的新平面
        \subitem （3）时序传播（Temporal Propagation）这部分是针对视频数据的优化，主要原理是认为时许接近的同名点具有相似性。
        \subitem （4）平面优化（Plane Refinement）。平面优化目标是通过优化视差平面$f_{p}$的参数而进一步减小像素p的聚合代价值$m(p,f_{p})$。PMS设置两个参数$\Delta^{max}_{z_{0}}$与$\Delta^{max}_{n}$。前者为点$P(x_{0},y_{0},z_{0})$的z坐标可变化范围，后者为平面法向量各分量的变化范围。在[-$\Delta^{max}_{z_{0}}$,$\Delta^{max}_{z_{0}}$]
        的范围内随机一个值$\Delta_{z_{0}}$加到$z_{0}$上得到$z^{'}_{0}=z_{0}+\Delta_{z_{0}}$。由此得到新的点。同样在法向量范围内随机一个法向量增量，得到新的法向量。新的点与新的法向量组成新的平面，进行代价聚合，如果新的平面代价低则更新像素p的平面。此过程也是迭代进行的，每次迭代就将变化范围缩小为原来的一半，直到迭代到特定范围内。
    \item \textbf{一致性检查与视差填充：} 在上述迭代过程结束后，将图像左右视图的视差进行比对，如果同名点视差在一定范围内则视为有效结果，否则结果无效。最后再检查图中视差的空洞对其进行填充。
    \end{enumerate}

\subsection{端到端深度学习MVS}
MVSNet 是端到端的有监督深度学习多视点深度图推断方法。MVSNet在传统的Plane-sweeping的基础上将代价体、特征点使用深度学习进行替换，从而构建了一个端到端的网络。
其主要流程及网络架构如图\ref{fig:MVSNetpipeline}所示。

\begin{figure}[H]
    \centering
    \includegraphics[height=0.3\textheight]{figures/MVSNet_pipeline.png}
    \caption{MVSNet流程}
    \label{fig:MVSNetpipeline}
\end{figure}

从图\ref{fig:MVSNetpipeline}中可以发现MVSNet主要包括了四个模块：

\begin{itemize}[leftmargin=4em]
    \item \textbf{特征提取模块：} 这部分主要是用CNN卷积神经网络将图像转换为高维特征图。所有图像共享一组CNN网络参数；
    \item \textbf{可微单应性投影：} 通过单应性变换将特征图投影到参考帧中，在这里论文使用的是可微单应性变换，使得通过深度学习框架自动求导得以实现，保持了整体框架的完整性，从而可以达到端到端的效果
    ，最后通过基于方差的代价计算构建代价体；
    \item \textbf{代价体正则化：}MVSNet 参考了UNet\cite[]{ronneberger2015u}，构建三维编码器与解码器。使用三维卷积对代价体进行正则化，使用softmax沿Channel维度进行运算得到归一化概率生成原始深度图；
    \item \textbf{深度图优化：} 通过引入参考图像对深度图进行优化。主要解决的是边界模糊的问题。
\end{itemize}

MVSNet作为基于深度学习的多视点深度估计和三维重建的开山之作，充分总结完善前人积淀，编码相机参数作为可微单应投影，并结合图像信息构建代价体，同时嫁接了 2D 图像
特征与3D损失回归的深度学习神经网络，为后续工作充分利用深度学习技术改进MVS问题的
求解打下了良好的基础。

另一方面，MVSNet作为有监督学习，对数据集有一定要求，获得带有真值的数据难度较大。同时端到端深度学习对各种环境的泛化能力也有待商榷。
本文将会在实验部分使用MVSNet在室内场景进行重建，观测其效果以论证这一点。

\subsection{本章小结}

本章按照传统三维重建的流程对其相关知识进行了介绍,并且分析了各种方法的优缺点。传统三维重建的精度较高,但在某些区域的效果不能让人满意。另一方面,本章介绍了端到端深度学习三维重建
的著名算法MVSNet,分析其原理与特点。本章为后续介绍本文方法提供了知识背景。